source("https://raw.githubusercontent.com/MatteoFasulo/Rare-Earth/main/R/util/coreFunctions.R")
loadPackages(c('sn','tidyverse','psych','RColorBrewer','stargazer','mclust','ContaminatedMixt',
'plotly','ggplot2','ggdendro','teigen','tclust','HDMD','caTools','clustvarsel',
'vscc','sparcl','pgmm','caret','glmnet','MLmetrics','DT'))
load("Z:\\DesktopC\\LUMSA\\2\\Data Mining\\Finite Mixture\\FiniteMixtureL31.RData")
rm(CO2data)
rm(NOdata)
rm(tonedata)
type <- wine$Type
rm(wine)Dati su 27 caratteristiche chimico/fisiche di tre diversi tipi di vino (Barolo, Grignolino, Barbera) dal Piemonte. Un set di dati con 178 osservazioni e 28 variabili (di cui la prima relativa alla tipologia di vino). Nell’ordine:
data(wines)
winesE’ stato possibile attraverso la ricerca originaria risalire all’anno di osservazione di ciascun vino. Di seguito vengono riportate le osservazioni dei tre diversi tipi di vino durante gli anni:
year <- as.numeric(substr(rownames(wines), 6, 7))
table(wines$wine, year) year
70 71 72 73 74 75 76 78 79
Barolo 0 19 0 20 20 0 0 0 0
Grignolino 9 9 7 9 16 9 12 0 0
Barbera 0 0 0 0 9 0 5 29 5
#wines[,'wine'] <- typeNotiamo subito che il Barbera è distribuito principalmente negli ultimi anni (76,78,79) mentre il Barolo nel 71, 73 e 74. Per quanto riguarda la percentuale delle singole classi:
wines %>%
count(wine = factor(wine)) %>%
mutate(pct = prop.table(n)) %>%
ggplot(aes(x = wine, y = pct, fill = wine, label = scales::percent(pct))) +
geom_col(position = 'dodge') +
geom_text(position = position_dodge(width = .9),
vjust = -0.5,
size = 3) +
scale_y_continuous(name = "Percentage")+
scale_x_discrete(name = "Wine Name")+
scale_fill_hue(l=40, c=35)+
theme(legend.position = "none")E’ chiaro che il Grignolino sia il più numeroso (39.9%) seguito dal Barolo (33.1%) e dal Barbera (27.0%).
Per rappresentare la dispersione dei dati abbiamo usato uno scatterplot leggermente differente dal solito. Sulla diagonale superiore si vede la distribuzione dei dati mentre sulla diagonale inferiore vi è la correlazione di Pearson tra le variabili.
# Correlation panel
my_cols <- c("#00AFBB", "#E7B800", "#FC4E07")
panel.cor <- function(x, y){
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- round(cor(x, y), digits=2)
txt <- paste0("R = ", r)
text(0.5, 0.5, txt, cex = 1)
}
# Customize upper panel
upper.panel<-function(x, y){
points(x,y, pch = 19, col = my_cols[wines$wine],cex=.5)
}
# Create the plots
pairs(wines[,2:8],
lower.panel = panel.cor,
upper.panel = upper.panel)Dalle prime 7 variabili, è possibile notare una correlazione di 0.69 tra acidity e malic e ovviamente una relazione inversamente proporzionale tra pH e acidity.
Per visualizzare le statistiche descrittive (media e deviazione standard) ci è sembrato opportuno dividerle in base alla classe di appartenenza:
printMeanAndSdByGroup <- function(variables,groupvariable)
{
variablenames <- c(names(groupvariable),names(as.data.frame(variables)))
groupvariable <- groupvariable[,1]
means <- aggregate(as.matrix(variables) ~ groupvariable, FUN = mean)
names(means) <- variablenames
print(paste("Means:"))
print(means)
sds <- aggregate(as.matrix(variables) ~ groupvariable, FUN = sd)
names(sds) <- variablenames
print(paste("Standard deviations:"))
print(sds)
}
printMeanAndSdByGroup(wines[2:28],wines[1])[1] "Means:"
[1] "Standard deviations:"
Alcune considerazioni:
Abbiamo calcolato la varianza within tra una feature e i tipi di vino:
calcWithinGroupsVariance(wines["flavanoids"],wines[1])[1] 0.2747075
Stesso discorso per la varianza between tra una feature e i vini:
calcBetweenGroupsVariance <- function(variable,groupvariable)
{
# find out how many values the group variable can take
groupvariable2 <- as.factor(groupvariable[[1]])
levels <- levels(groupvariable2)
numlevels <- length(levels)
# calculate the overall grand mean:
grandmean <- sapply(variable,mean)
# get the mean and standard deviation for each group:
numtotal <- 0
denomtotal <- 0
for (i in 1:numlevels)
{
leveli <- levels[i]
levelidata <- variable[groupvariable==leveli,]
levelilength <- length(levelidata)
# get the mean and standard deviation for group i:
meani <- mean(levelidata)
sdi <- sd(levelidata)
numi <- levelilength * ((meani - grandmean)^2)
denomi <- levelilength
numtotal <- numtotal + numi
denomtotal <- denomtotal + denomi
}
# calculate the between-groups variance
Vb <- numtotal / (numlevels - 1)
Vb <- Vb[[1]]
return(Vb)
}
calcBetweenGroupsVariance(wines["flavanoids"],wines[1])[1] 64.2612
Per vedere quali variabili hanno la maggiore separazione, (rapporto tra Varianza Between e Varianza Within) abbiamo scritto una funzione apposita per calcolarne il valore per ogni feature.
Per alcune analisi successive, abbiamo provato a cambiare il task della ricerca tentando di utilizzare un modello come classificatore e di misurarne le prestazioni di classificazione. A tale scopo, abbiamo suddiviso il dataset originario in due sottogruppi:
require(caTools)
sample = sample.split(wines[,1], SplitRatio = .50)
train = subset(wines, sample == TRUE)
trainTestNames <- train$wine
print(paste("Train Obs:",nrow(train)))[1] "Train Obs: 90"
#train$wine <- as.numeric(train$wine)
test = subset(wines, sample == FALSE)
wineTestNames <- test$wine
print(paste("Test Obs:",nrow(test)))[1] "Test Obs: 88"
#test$wine <- as.numeric(test$wine)Per il Clustering abbiamo deciso di applicare:
dissMatrix <- pairwise.mahalanobis(wines[,-1],
grouping = c(1:nrow(wines)),
cov = cov(wines[,-1]))$distance
dissMatrix <- sqrt(dissMatrix)
dissMatrix <- as.dist(dissMatrix)
combDist <- function(distance, methods, df, dt, dissMatrix) {
c <- 0
results <- list()
for (i in 1:length(distance)){
ifelse(distance[i] == "minkowski",
dist <- dist(df, method = distance[i], p = 4),
ifelse(distance[i] == "mahalanobis",
dist <- dissMatrix,
dist <- dist(df, method = distance[i])))
for (j in 1:length(methods)){
dendo = hclust(dist, method = methods[j])
dendo.x = ggdendrogram(dendo, rotate = F, size = 2, leaf_labels = T, labels = F) +
ggtitle(paste(distance[i],' ',methods[j],sep=''))
for(elem in 2:4){
cluster = cutree(dendo, k=elem)
c <- c + 1
results[[c]] <- list(distance = distance[i],
method = methods[j],
groups = elem,
table = table(dt,cluster),
dendo = dendo.x,
AdjustedRandIndex = adjustedRandIndex(dt,cluster),
cluster = cluster)
}
}
}
return(results)
}
results <- combDist(c("euclidean", "manhattan", "minkowski","mahalanobis"),
c("single", "complete", "average", "ward.D"), scale(wines[,-1]), wines[,1], dissMatrix)
optimal <- function(results){
best_randIndex.eu = 0
best_randIndex.ma = 0
best_randIndex.mi = 0
best_randIndex.maha = 0
best_model.eu = integer()
best_model.ma = integer()
best_model.mi = integer()
best_model.maha = integer()
for (i in 1:length(results)){
current_randIndex = results[[i]]$AdjustedRandIndex
if (results[[i]]$distance == "euclidean"){
if (current_randIndex > best_randIndex.eu) {
best_randIndex.eu = current_randIndex
best_model.eu = i
}
}
else if (results[[i]]$distance == "manhattan"){
if (current_randIndex > best_randIndex.ma) {
best_randIndex.ma = current_randIndex
best_model.ma = i
}
}
else if (results[[i]]$distance == "minkowski"){
if (current_randIndex > best_randIndex.mi) {
best_randIndex.mi = current_randIndex
best_model.mi = i
}
}
else if (results[[i]]$distance == "mahalanobis"){
if (current_randIndex > best_randIndex.maha) {
best_randIndex.maha = current_randIndex
best_model.maha = i
}
}
}
#print(list(euclidean = list(model.number = best_model.eu,
# cluster = results[[best_model.eu]]$groups,
# AdjustedRandIndex = best_randIndex.eu),
# manhattan = list(model.number = best_model.ma,
# cluster = results[[best_model.ma]]$groups,
# AdjustedRandIndex = best_randIndex.ma),
# minkowski = list(model.number = best_model.mi,
# cluster = results[[best_model.mi]]$groups,
# AdjustedRandIndex = best_randIndex.mi),
# mahalanobis=list(model.number = best_model.maha,
# cluster = results[[best_model.maha]]$groups,
# AdjustedRandIndex = best_randIndex.maha))
# )
return(list(euclidean = results[[best_model.eu]],
manhattan = results[[best_model.ma]],
minkowski = results[[best_model.mi]],
mahalanobis=results[[best_model.maha]]))
}
best_dist_model = optimal(results)print(paste("AdjustedRandIndex:",round(best_dist_model$euclidean$AdjustedRandIndex,3)))[1] "AdjustedRandIndex: 0.769"
print(paste("AdjustedRandIndex:",round(best_dist_model$manhattan$AdjustedRandIndex,3)))[1] "AdjustedRandIndex: 0.946"
ggplotly(best_dist_model$minkowski$dendo)print(best_dist_model$minkowski$table) cluster
dt 1 2 3 4
Barolo 55 2 2 0
Grignolino 8 1 59 3
Barbera 2 0 0 46
print(paste("AdjustedRandIndex:",round(best_dist_model$minkowski$AdjustedRandIndex,3)))[1] "AdjustedRandIndex: 0.73"
print(paste("AdjustedRandIndex:",round(best_dist_model$mahalanobis$AdjustedRandIndex,4)))[1] "AdjustedRandIndex: 0.27"
require(cluster)
k.means.3 <- kmeans(scale(wines[,-1]),centers=3,nstart = 50, iter.max = 100)
table(wines[,1], k.means.3$cluster)
1 2 3
Barolo 1 0 58
Grignolino 68 1 2
Barbera 0 48 0
print(paste("AdjustedRandIndex:",round(adjustedRandIndex(k.means.3$cluster, wines[,1]),3)))[1] "AdjustedRandIndex: 0.93"
require(cluster)
k.means.4 <- kmeans(scale(wines[,-1]),centers=4,nstart = 50, iter.max = 100)
table(wines[,1], k.means.4$cluster)
1 2 3 4
Barolo 0 0 59 0
Grignolino 38 0 5 28
Barbera 0 47 0 1
print(paste("AdjustedRandIndex:",round(adjustedRandIndex(k.means.4$cluster, wines[,1]),3)))[1] "AdjustedRandIndex: 0.736"
require(cluster)
PAM.3 <- pam(wines[,-1], k=3,
metric = "euclidean",
nstart = 50,
stand = TRUE)
table(wines[,1], PAM.3$clustering)
1 2 3
Barolo 50 9 0
Grignolino 3 67 1
Barbera 1 1 46
print(paste("AdjustedRandIndex:",round(adjustedRandIndex(PAM.3$clustering, wines[,1]),3)))[1] "AdjustedRandIndex: 0.754"
require(cluster)
PAM.4 <- pam(wines[,-1], k=4,
metric = "euclidean",
nstart = 50,
stand = TRUE)
table(wines[,1], PAM.4$clustering)
1 2 3 4
Barolo 28 27 4 0
Grignolino 3 0 67 1
Barbera 1 0 1 46
print(paste("AdjustedRandIndex:",round(adjustedRandIndex(PAM.4$clustering, wines[,1]),3)))[1] "AdjustedRandIndex: 0.728"
Novacula Occami: frustra fit per plura quod potest fieri per pauciora (Il rasoio di Occam: è futile fare con più mezzi ciò che si può fare con meno). Tale principio metodologico è ritenuto alla base del pensiero scientifico moderno.
1 2 3 4
Barolo 47 10 2 0
Grignolino 6 3 61 1
Barbera 0 0 1 47
[1] "AdjustedRandIndex: 0.734"
1 2 3 4
Barolo 47 10 2 0
Grignolino 6 3 61 1
Barbera 0 0 1 47
[1] "AdjustedRandIndex: 0.734"
table(wines[,1], vscc.mclust$initialrun$classification) #Clustering results on full data set
1 2 3
Barolo 58 1 0
Grignolino 1 70 0
Barbera 0 2 46
print(paste("AdjustedRandIndex:",round(adjustedRandIndex(vscc.mclust$initialrun$classification, wines[,1]),3)))[1] "AdjustedRandIndex: 0.931"
table(wines[,1], vscc.mclust$bestmodel$classification) #Clustering results on reduced data set
1 2 3
Barolo 59 0 0
Grignolino 3 66 2
Barbera 0 0 48
print(paste("AdjustedRandIndex:",round(adjustedRandIndex(vscc.mclust$bestmodel$classification, wines[,1]),3)))[1] "AdjustedRandIndex: 0.913"
mixt.wines <- Mclust(wines[,-1],G=3:8)
cn.wines.mixt <- CNmixt(wines[,-1], G = 3, initialization = "mixt", seed = 1234, parallel = F, verbose = F)
cn.wines.kmeans <- CNmixt(wines[,-1], G = 3, initialization = "kmeans", seed = 1234, parallel = F, verbose = F)
cn.wines.rpost <- CNmixt(wines[,-1], G = 3, initialization = "random.post", seed = 1234, parallel = F, verbose = F)
cn.wines.rclass <- CNmixt(wines[,-1], G = 3, initialization = "random.clas", seed = 1234, parallel = F, verbose = F)
teigen.kmeans <- teigen(wines[,-1], Gs=3:4, init = 'kmeans', scale = T)
adjustedRandIndex(wines[,1],teigen.kmeans$classification)
teigen.classifier.kmeans <- teigen(train[,-1], Gs=3:4, init = 'kmeans', scale = T, known = train[,1])
predict(teigen.classifier.kmeans,test[,-1])
adjustedRandIndex(test[,1],predict(teigen.classifier.kmeans, test[,-1])$classification)
teigen.classifier.uniform <- teigen(train[,-1], Gs=3:4, init = 'uniform', scale = T, known = train[,1])
predict(teigen.classifier.uniform,test[,-1])
adjustedRandIndex(test[,1],predict(teigen.classifier.uniform, test[,-1])$classification)
teigen.hard <- teigen(wines[,-1], Gs=3:4, init = 'hard', scale = T)
adjustedRandIndex(wines[,1],teigen.hard$classification)
teigen.soft <- teigen(wines[,-1], Gs=3:4, init = 'soft', scale = T)
adjustedRandIndex(wines[,1],teigen.soft$classification)Per la nostra analisi abbiamo voluto verificare l’efficienza di tre noti modelli di regolarizzazione attraverso il pacchetto caret:
lambda <- 10^seq(0, -2, length = 250)Il modello Ridge riduce i coefficienti, in modo che le variabili, con un contributo minore al risultato, abbiano i loro coefficienti vicini allo zero. Invece di forzarli a essere esattamente zero (come nel Lasso), li penalizziamo con un termine chiamato norma L2 costringendoli così a essere piccoli. In questo modo diminuiamo la complessità del modello senza eliminare nessuna variabile attraverso una costante chiamata lambda (\(\lambda\)) di penalizzazione: \[ L_{ridge}(\hat{\beta}) = \sum_{i = 1}^{n}{(y_i - x_i\hat{\beta})^2} + \lambda\sum_{k = 1}^{K}{\hat{\beta}_k^2} \]
# Build the model
set.seed(123)
ridge <- caret::train(
x = train[,-1],
y = factor(train[,1]),
method = "glmnet",
trControl = trainControl("cv", number = 10, classProbs = TRUE, summaryFunction = multiClassSummary),
tuneGrid = expand.grid(alpha = 0, lambda=lambda),
metric="Accuracy")
# Model coefficients
coef(ridge$finalModel, ridge$bestTune$lambda)
# Make predictions
predictions.ridge <- ridge %>% predict(test)
# Model prediction performance
tibble(
trueValue = wineTestNames,
predictedValue = predictions.ridge)Il ridge è composto dalla somma dei residui quadrati più una penalità, definita dalla lettera Lambda, che è moltiplicata per la somma dei coefficienti quadrati \(\beta\). Quando \(\lambda = 0\), il termine di penalità non ha alcun effetto e il ridge produrrà i coefficienti minimi quadrati classici. Tuttavia, quando \(\lambda\) aumenta all’infinito, l’impatto della penalità aumenta e i coefficienti si avvicinano allo zero. Il ridge è particolarmente indicato quando si hanno molti dati multivariati con numero di feature maggiore del numero di osservazioni. Lo svantaggio, però, è che includerà tutti le feature nel modello finale, a differenza dei metodi di feature selection, che generalmente selezioneranno un insieme ridotto di variabili tra quelle disponibili.
caret::confusionMatrix(predictions.ridge, test$wine)Il Least Absolute Shrinkage and Selection Operator (LASSO) riduce i coefficienti verso lo zero penalizzando il modello con un termine di penalità chiamato norma L1, che è la somma dei coefficienti in valore assoluto: \[ L_{lasso}(\hat{\beta}) = \sum_{i = 1}^{n}{(y_i - x_i\hat{\beta})^2} + \lambda\sum_{k = 1}^{K}{|\hat{\beta}_k|} \]
# Build the model
set.seed(123)
lasso <- caret::train(
x = train[,-1],
y = factor(train[,1]),
method = "glmnet",
trControl = trainControl("cv", number = 10, classProbs = TRUE, summaryFunction = multiClassSummary),
tuneGrid = expand.grid(alpha = 1, lambda=lambda),
metric="Accuracy")
# Model coefficients
coef(lasso$finalModel, lasso$bestTune$lambda)
# Make predictions
predictions.lasso <- lasso %>% predict(test)
# Model prediction performance
tibble(
trueValue = wineTestNames,
predictedValue = predictions.lasso)In questo caso la penalità ha l’effetto di forzare alcune delle stime dei coefficienti, con un contributo minore al modello, a essere esattamente uguale a zero. Il lasso, quindi, può anche essere visto come un’alternativa ai metodi di feature selection per eseguire la selezione delle variabili al fine di ridurre la complessità del modello.
Come nel ridge, è fondamentale selezionare un buon valore di \(\lambda\).
Quando lambda è piccolo, il risultato è molto vicino alla stima dei minimi quadrati. All’aumentare di lambda, si verifica una contrazione in modo da poter eliminare le variabili che sono a zero.
caret::confusionMatrix(predictions.lasso, test$wine)Elastic Net combina le proprietà di Ridge e Lasso penalizzando il modello usando sia la norma L2 che la norma L1: \[ L_{ElasticNet}(\hat{\beta}) = \frac{\sum_{i = 1}^{n}{(y_i - x_i\hat{\beta})^2}}{2n} + \lambda(\frac{1-\alpha}{2}\sum_{k = 1}^{K}{\hat{\beta}_k^2} + \alpha\sum_{k = 1}^{K}{|\hat{\beta}_k}|) \]
set.seed(123)
elastic <- caret::train(
x = train[,-1],
y = factor(train[,1]),
method = "glmnet",
trControl = trainControl("cv", number = 10, classProbs = TRUE, summaryFunction = multiClassSummary),
metric="Accuracy")
# Model coefficients
coef(elastic$finalModel, elastic$bestTune$lambda)
# Make predictions
predictions.enet <- elastic %>% predict(test)
# Model prediction performance
tibble(
trueValue = wineTestNames,
predictedValue = predictions.enet)Oltre a impostare e scegliere un valore lambda, l’elastic net ci consente anche di ottimizzare il parametro alfa dove \(\alpha = 0\) corrisponde a ridge e \(\alpha = 1\) al lasso.
Pertanto possiamo scegliere un valore \(\alpha\) compreso tra 0 e 1 per ottimizzare l’elastic net. Se tale valore è incluso in questo intervallo, si avrà una riduzione con alcuni portati a \(0\).
caret::confusionMatrix(predictions.enet, test$wine)models <- list(ridge = caret::confusionMatrix(predictions.ridge, test$wine),
lasso = caret::confusionMatrix(predictions.lasso,test$wine),
elastic = caret::confusionMatrix(predictions.enet,test$wine))
#models %>% summary(metric = "Accuracy")# Il rilassamento dato da (p-q)^2 > p+q è verificato per 1 <= q <= 20
pg.random.ccc <- pgmmEM(selectedWines[,-1], rG=3:4, rq=1:5, seed = 1234, relax = T)
adjustedRandIndex(wines[,1],pg.random.ccc$map)
pg.random <- pgmmEM(wines[,-1], rG=3:4, rq=20:25,icl=F,zstart = 1,cccStart = F,loop = 50, seed = 1234, relax = T)
adjustedRandIndex(wines[,1],pg.random.ccc$map)
pg2 <- pgmmEM(wines[,-1], rG=3:4, rq=20:25,icl=T,zstart = 1,cccStart = T,loop = 50, seed = 1234, relax = T)
adjustedRandIndex(wines[,1],pg2$map)
pg2 <- pgmmEM(wines[,-1], rG=3:4, rq=20:25,icl=T,zstart = 2, seed = 1234, relax = T)
adjustedRandIndex(wines[,1],pg2$map)
pg2 <- pgmmEM(wines[,-1], rG=3:4, rq=20:25,icl=T,zstart = 3, seed = 1234, relax = T, zlist=)
adjustedRandIndex(wines[,1],pg2$map)
agree(cn.wines.mixt, givgroup = wines[,1])
agree(cn.wines.kmeans, givgroup = wines[,1])
agree(cn.wines.rpost, givgroup = wines[,1])
agree(cn.wines.rclass, givgroup = wines[,1])
plot(cn.wines.mixt, contours = TRUE)
plot(cn.wines.kmeans, contours = TRUE)
plot(cn.wines.rpost, contours = TRUE)
plot(cn.wines.rclass, contours = TRUE)if (class(models[[i]])==‘kmeans’){ ariTrain[i] <- adjustedRandIndex(train,models[[i]]\(cluster) ariTest[i] <- adjustedRandIndex(test,models[[i]]\)cluster) bic[i] <- NA icl[i] <- NA G[i] <- length(unique(models[[i]]$cluster))
summaryOfModels <- function(train, test, models){
nModel <- length(models)
ariTrain <- NULL
ariTest <- NULL
bic <- NULL
icl <- NULL
G <- NULL
for(i in 1:nModel){
if (class(models[[i]])=='Mclust'){
ariTrain[i] <- adjustedRandIndex(train,models[[i]]$classification)
ariTest[i] <- adjustedRandIndex(test,predict(models[[i]],test)$classification)
bic[i] <- models[[i]]$bic
icl[i] <- models[[i]]$icl
G[i] <- models[[i]]$G
} else if (class(models[[i]])=='ContaminatedMixt'){
bestModel <- getBestModel(models[[i]])$models[[1]]
ariTrain[i] <- adjustedRandIndex(train,bestModel$group)
ariTest[i] <- adjustedRandIndex(test,predict(bestModel,))
bic[i] <- whichBestModel(models[[i]],train)$BIC
icl[i] <- whichBestModel(models[[i]],train)$ICL
G[i] <- whichBestModel(models[[i]],train)$G
} else if (class(models[[i]])=='teigen'){
ariTrain[i] <- adjustedRandIndex(train,models[[i]]$iclresults$classification)
ariTest[i] <- adjustedRandIndex(test,models[[i]]$iclresults$classification)
bic[i] <- models[[i]]$bic
icl[i] <- models[[i]]$iclresults$icl
G[i] <- models[[i]]$G
} else if (class(models[[i]])=='pgmm'){
ari[i] <- adjustedRandIndex(df,models[[i]]$map)
G[i] <- models[[i]]$g
ifelse(is.null(models[[i]]$bic[1])==TRUE,
bic[i] <- NA,
bic[i] <- as.double(models[[i]]$bic[1]))
ifelse(is.null(models[[i]]$icl[1])==TRUE,
icl[i] <- NA,
icl[i] <- as.double(models[[i]]$icl[1]))
} else if (class(models[[i]])=='tkmeans'){
ari[i] <- adjustedRandIndex(df,models[[i]]$cluster)
G[i] <- models[[i]]$k
bic[i] <- NA
icl[i] <- NA
}
}
outputDF <- data.frame('AdjustedRandIndex' = ari,
'BIC' = bic,
'ICL' = icl,
'G' = as.integer(G),
row.names = c('KMeans','Mclust','ContaminatedMixt','TEigen','PGMM','TClust'),
stringsAsFactors = F)
return(outputDF)
}
test <- summaryOfModels(wines[,1],list(k.means,mixt.wines,mixt.cn.wines,teigen_wine,pippo,trimmedOne))
test